Deep Learning Analysis of Defect and Phase Evolution During Electron Beam Induced Transformations in WS$_2$

Authors: Artem Maksov, Maxim Ziatdinov*

*Correspondence to: ziatdinovmax@gmail.com

Date: 03/09/2018

In this notebook, we demonstrate a complete workflow for extraction and classification of defects in a STEM image stack (STEM "movie").

Table of Contents

  1. We use a single frame from a movie with FFT subtraction method to create a training dataset.
  2. We train a fully convolutional neural network (without "dense" layers) to find defects that break lattice periodicity.
  3. We apply the network to the whole movie and extract subwindows containing individual defects.
  4. We train a Gaussian Mixture Model to classify the defects and then we explore their structure and evolution in time.

Requirenments:

  1. Python 3 (matplolib, numpy, scipy, sklearn, opencv, h5py)
  2. Latest version of pydm3reader: https://microscopies.med.univ-tours.fr/?page_id=260
  3. Keras + tensorflow backend

0. Imports

In [1]:
%matplotlib nbagg
#%matplotlib inline

#I/O
import os
import dm3_lib as dm3
import h5py
from sklearn.externals import joblib

#graphics
import matplotlib
import matplotlib.pylab as plt
import matplotlib.image as mpimg
import matplotlib.patches as patches
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Line3DCollection
import matplotlib.gridspec as gridspec
import pylab as P


#processing 
import numpy as np
from scipy import fftpack
from scipy import ndimage

from skimage.feature import blob_log
from sklearn import cluster, mixture
from scipy.cluster.hierarchy import dendrogram, linkage
from sklearn.decomposition import PCA
from collections import OrderedDict
from scipy.spatial import cKDTree
In [2]:
import cv2
#neural networks
from keras.layers import Input, Convolution2D, MaxPooling2D, UpSampling2D
from keras.layers import Dropout, Activation, Reshape
from keras.models import Model, load_model
from keras.callbacks import TensorBoard
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import OneHotEncoder
Using TensorFlow backend.

1. Initial image processing

1.1. Reading in the file

Define path to the file:

In [3]:
path = './'
filename = 'stack4.dm3'

Read in the image data:

In [4]:
dm3f = dm3.DM3(path + filename)
imgdata = dm3f.imagedata

We need to rescale the image according to gain, and normalize it:

In [5]:
for t in dm3f.tags:
    if 'PMTDF' in t:
        scale = float(dm3f.tags[t])/64000
        
imgdata = imgdata/scale
In [6]:
imgdata = imgdata - np.amin(imgdata)
imgdata = imgdata/np.amax(imgdata)

Create directory for storing outputs:

In [9]:
directory = './Figures/' + filename[:-4] + '/'
if not os.path.exists(directory):
    os.makedirs(directory)

1.2 FFT subtraction

We will explore the first frame of the movie to estimate parameters for FFT subtraction

In [10]:
imgsrc = imgdata[0]

1.2.1 Plotting real space image:

In [63]:
#plot real space image
fig100 = plt.figure(100, figsize=(5,5))
ax = fig100.add_axes([0, 0, 1, 1])
ax.imshow(imgsrc, cmap = 'gray')
ax.axis('off')

#save the figure
fig100.savefig(directory + filename + '_realspace.png')

1.2.2 FFT subtraction

We perform FFT subtraction in a similar manner as to described in (reference)

We need to create a mask and apply it as a high-pass filter in the frequency domain. Currently, this depends on 1 human-defined parameter: 1/ratio of radius of the mask to the image size.

In [64]:
# 1/ratio of the radius of the mask to the original image in pixels
maskratio = 10
In [65]:
def FFTmask(imgsrc, maskratio=10):
    """Takes a square real space image and filter out a disk with radius equal to:
    1/maskratio * image size.
    Retruns FFT transform of the image and the filtered FFT transform
    """
    # Take the fourier transform of the image.
    F1 = fftpack.fft2((imgsrc)) 
    # Now shift so that low spatial frequencies are in the center.
    F2 = (fftpack.fftshift((F1)))
    # copy the array and zero out the center
    F3 = F2.copy()
    l = int(imgsrc.shape[0]/maskratio)
    m = int(imgsrc.shape[0]/2)
    y,x = np.ogrid[1: 2*l + 1, 1:2*l + 1]
    mask = (x - l)*(x - l) + (y - l)*(y - l) <= l*l
    F3[m-l:m+l, m-l:m+l] = F3[m-l:m+l, m-l:m+l] * (1 - mask)
    
    return F2, F3
In [66]:
F2, F3 = FFTmask(imgsrc, maskratio)

We can plot both FFT and filtered FFT to see how it performs:

In [67]:
#plot 
fig101 = plt.figure(101)
ax1 = fig101.add_subplot(121)
ax1.imshow(np.abs(np.log(F2)), cmap='hot')
ax1.axis('off')
ax2 = fig101.add_subplot(122)
masked_array = np.ma.array (np.abs(np.log(F3)), mask=np.isnan(np.abs(np.log(F3))))
cm = matplotlib.cm.hot
cm.set_bad('black',1.)
ax2.imshow(masked_array, cm)
ax2.axis('off')

fig101.savefig(directory + 'FFTfilter.png')
c:\anaconda3\envs\tensorflow\lib\site-packages\ipykernel_launcher.py:7: RuntimeWarning: divide by zero encountered in log
  import sys

Now we reconstruct the real space image from filtered FFT image and subtract it from the original to identify locations with broken symmetry

In [20]:
def FFTsub(imgsrc, F3):
    """Takes real space image and filtred FFT.
    Reconstructs real space image and subtracts it from the original.
    Returns normalized image. 
    """
    reconstruction = np.real(fftpack.ifft2(fftpack.ifftshift(F3)))
    diff = np.abs(imgsrc - reconstruction)
    
    #normalization
    diff = diff - np.amin(diff)
    diff = diff/np.amax(diff)
    
    return diff
In [21]:
diff = FFTsub(imgsrc, F3)
In [22]:
fig102 = plt.figure(102, figsize=(5,5))
ax = fig102.add_axes([0, 0, 1, 1])
ax.imshow(diff, cmap='gray')
ax.axis('off')
fig102.savefig(directory + 'diff.png')

1.2.3 Creating a defect location map

We need to thresold the resulting difference image to create a map for the defect locations, which will serve as labels for pixels in our training set

We need to define 2 thresholds: low for low intensity defects and high for high-intensity. Since we normalized the difference image they are both in [0;1) range

We can use histogram of intensities to inform our choice:

In [23]:
diffall = diff.reshape(diff.shape[0] * diff.shape[1])
In [24]:
fig103 = plt.figure(103)
ax103 = fig103.add_axes([0, 0, 1, 1])
ax103.hist(diffall, bins='auto');

Or we can inform our choice by computing statistics for the intensity distributions, but those should be taken as guidance and confirmed by visual inspection as they can vary from image to image

In [25]:
print('Defining cutoff as 2 std, low threshold:', np.mean(diffall) - 2* np.std(diffall), 
      ' high thresold: ', np.mean(diffall) + 2* np.std(diffall))
print('Defining cutoff as 3 std, low threshold:', np.mean(diffall) - 3* np.std(diffall), 
      ' high thresold: ', np.mean(diffall) + 3* np.std(diffall))
Defining cutoff as 2 std, low threshold: 0.30713032186  high thresold:  0.631990864873
Defining cutoff as 3 std, low threshold: 0.225915186107  high thresold:  0.713206000626

Define the thresholds:

In [27]:
def threshImg(diff, threshL=0.25, threshH=0.75):
    """Takes in difference image, low and high thresold values, and outputs a map of all defects.
    """
    
    threshIL = diff < threshL  
    threshIH = diff > threshH
    threshI = threshIL + threshIH
    
    return threshI
In [31]:
threshL = 0.25
threshH = 0.75

Plot the results:

In [32]:
fig104 = plt.figure(104, figsize=(5,5))
ax = fig104.add_axes([0, 0, 1, 1])

threshI = threshImg(diff, threshL, threshH)

ax.imshow(imgsrc)
ax.imshow(threshI, cmap='gray', alpha=0.3)
ax.axis('off')
Out[32]:
(-0.5, 511.5, 511.5, -0.5)
In [33]:
labelIm = threshI.astype(int)

1.2.4. Generating training set

We can generate training set by taking patches of the image and corresponding labels

We define helper functions:

In [34]:
#Generate x, y positions for sliding windows
def GenerateXYPos(window_size, window_step, image_width):
    """Takes the window size, step, and total image width 
    and generates all xy pairs for sliding window"""   
    xpos_vec = np.arange(0,image_width-window_size,window_step)
    ypos_vec = np.arange(0,image_width-window_size,window_step)
        
    num_steps = len(xpos_vec)    
            
    xpos_mat = np.tile(xpos_vec, num_steps)
    ypos_mat = np.repeat(ypos_vec, num_steps)
    pos_mat = np.column_stack((xpos_mat, ypos_mat))
            
    return pos_mat
In [35]:
def MakeWindow(imgsrc, xpos, ypos, window_size):
    """returns window of given size taken at position on image"""
    imgsrc = imgsrc[xpos:xpos+window_size, ypos:ypos+window_size]
    return imgsrc
In [36]:
def imgen(raw_image, pos_mat, window_size):
    """Returns all windows from image for given positions array"""
    
    immat = np.zeros(shape = (len(pos_mat), window_size, window_size))
    
    for i in np.arange(0,len(pos_mat)):
        img_window = MakeWindow(raw_image, pos_mat[i,0], pos_mat[i,1], window_size) 
        immat[i,:,:,] = img_window
    
    return immat
In [37]:
def image_preprocessing(image_data, norm=0):
    """Reshapes data and optionally normalizes it"""
    image_data = image_data.reshape(image_data.shape[0], image_data.shape[1], image_data.shape[2], 1)
    image_data = image_data.astype('float32')
    if norm != 0:
        image_data = (image_data - np.amin(image_data))/(np.amax(image_data) - np.amin(image_data))
    return image_data
In [41]:
def label_preprocessing(image_data, nb_classes):
    """Returns labels / ground truth for images"""
    
    label4D = np.empty((0, image_data.shape[1], image_data.shape[2], nb_classes))
    for idx in range(image_data.shape[0]):
        img = image_data[idx,:,:]
        n, m = img.shape
        img = np.array(OneHotEncoder(n_values=nb_classes).fit_transform(img.reshape(-1,1)).todense())
        img = img.reshape(n, m, nb_classes)
        label4D = np.append(label4D, [img], axis = 0)    
    return label4D
In [ ]:
def label_preprocessing2(image_data):
    """we can simplify this in case of only two classes"""
    
    label1 = image_data.reshape(image_data.shape[0], image_data.shape[1], image_data.shape[2], 1)
    label2 = -label1 + 1
    label4D = np.concatenate((label2, label1), axis = 3)
    return label4D

We can control the number of patches by changing window size and window step

In [42]:
window_size = 256
window_step = 8
pos_mat = GenerateXYPos(window_size, window_step, imgsrc.shape[0])
print("On image size: ", imgsrc.shape[0], ' by ',  imgsrc.shape[1], ' pixels: ', len(pos_mat), "patches will be generated")
On image size:  512  by  512  pixels:  1024 patches will be generated
In [43]:
immat = imgen(imgsrc, pos_mat, window_size)
labelmat = imgen(labelIm, pos_mat, window_size)
In [44]:
immat2 =  image_preprocessing(immat, 1)
labelmat2 =  label_preprocessing(labelmat, nb_classes = 2)

2. Training the Neural Network model

First, we need to define network architecture:

In [46]:
def model_defect(input_img, nb_classes = 2):
    """Creates a Deep Learning model for defect identification"""
    
    x = Convolution2D(20, (3, 3), activation='relu', padding='same')(input_img)
    x = MaxPooling2D((2, 2), padding='same')(x)
    x = Convolution2D(40, (3, 3), activation='relu', padding='same')(x)
    x = MaxPooling2D((2, 2))(x)
        
    x = Convolution2D(40, (3, 3), activation='relu', padding='same')(x)
    x = UpSampling2D((2, 2))(x)
    x = Convolution2D(20, (3, 3), activation='relu', padding='same')(x)
    x = UpSampling2D((2, 2))(x)
    
    x = Convolution2D(nb_classes, (3, 3), activation = 'linear', padding='same')(x)
    x = Convolution2D(nb_classes, (1, 1), activation = 'linear', padding='same')(x)

    output = Activation('softmax')(x)
    
    return Model(input_img, output)

We can split data into training and validation set to improve training:

In [52]:
nb_classes = 2
target_size = (256, 256)

X_train, X_test, y_train, y_test = train_test_split(immat2, labelmat2, test_size = 0.2, random_state = 42)

Finally, we train the model:

In [ ]:
input_img = Input(shape=(target_size[0], target_size[1], 1)) 

model = model_defect(input_img)
model.compile(optimizer='adam', loss='categorical_crossentropy', metrics = ['accuracy'])
callback_tb = TensorBoard(log_dir='/tmp/AtomGen', histogram_freq=0,
                          write_graph=True, write_images=False)

model.fit(X_train, y_train, epochs=50, batch_size=32, 
          validation_data=(X_test, y_test), shuffle=True, verbose = 1)

And save it:

In [216]:
model.save("./DefectModel3.h5")
model.save_weights("./DefectModel_weights3.h5")
print('Saved model and weights to disk.\n')
Saved model and weights to disk.

3. Processing the movie

We can apply the trained model to the rest of the movie to identify defects.

3.1 Processing single frame

In [8]:
nframes = imgdata.shape[0]
In [58]:
def splitImage(img, target_size):
    """Splits image into patches with given size"""
    
    xs = img.shape[0]
    ys = img.shape[1]
    nxp = int(xs/target_size[0])
    nyp = int(xs/target_size[1])
    
    impatchmat = np.zeros(shape=(int(nxp*nyp), target_size[0], target_size[1], 1))
    
    count = 0
    for i in range(nxp):
        for i2 in range(nyp):
            xstart = target_size[0] * i
            ystart = target_size[1] * i2
            xend = target_size[0] * (i + 1)
            yend = target_size[1] * (i2 + 1)

            impatchmat[count, :, :, 0] = img[xstart:xend, ystart:yend]
            count = count + 1
            
    return impatchmat
In [80]:
def predictDefects(img, model, target_size, nb_classes=2):
    """Uses given DL model to generate prediciton maps on the image"""
    
    xs = img.shape[0]
    ys = img.shape[1]
    nxp = int(xs/target_size[0])
    nyp = int(xs/target_size[1])
    classpred = np.zeros(shape=(nb_classes,xs, ys))

    impatchmat = splitImage(img, target_size)
    res = model.predict(impatchmat)
    
    count = 0
    for i in range(nxp):
        for i2 in range(nyp):
            xstart = target_size[0] * i
            ystart = target_size[1] * i2
            xend = target_size[0] * (i + 1)
            yend = target_size[1] * (i2 + 1)
            
            for i3 in range(nb_classes):
                classpred[i3, xstart:xend, ystart:yend] = res[count, :, :, i3]
                
            count = count + 1
            
    return classpred

Here is just a test for a single frame:

In [129]:
fnum = 80
classpred = predictDefects(imgdata[fnum], model, target_size, nb_classes)
In [130]:
fig300 = plt.figure(300)
ax1 = fig300.add_subplot(121)
ax1.imshow(imgdata[fnum])
ax1.axis('off')
ax2 = fig300.add_subplot(122)
ax2.imshow(classpred[1])
ax2.axis('off')
Out[130]:
(-0.5, 511.5, 511.5, -0.5)
In [131]:
fig301 = plt.figure(301, figsize=(5,5))
ax = fig301.add_axes([0, 0, 1, 1])

ax.imshow(imgdata[fnum])
ax.imshow(classpred[1] > 0.99, cmap='hot', alpha=0.5)
ax.axis('off')
Out[131]:
(-0.5, 511.5, 511.5, -0.5)

3. 2 Defect location identification

To classify the defects, we need to extract image patches containing them. This means we need to filter out defects with large areas, which can be observed in the second half of the movie.

In [68]:
softmax_threhold = 0.5
In [140]:
def extractDefects(img, classpred, softmax_threhold = 0.5, bbox=16):
    
    defim = np.ndarray(shape=(0, bbox*2, bbox*2))
    defcoord = np.ndarray(shape=(0, 2))
    defcount = 0
    
    _, thresh = cv2.threshold(classpred[1], softmax_threhold, 1, cv2.THRESH_BINARY)
    
    s = [[1,1,1],[1,1,1],[1,1,1]]
    labeled, nr_objects = ndimage.label(thresh, structure=s) 
    loc = ndimage.find_objects(labeled)
    cc = ndimage.measurements.center_of_mass(labeled, labeled, range(nr_objects + 1))
    sizes = ndimage.sum(thresh,labeled,range(1, nr_objects +1)) 

    #filter found points
    cc2 = cc[1:]
    t = zip(cc2, sizes)
    ccc = [k[0] for k in t if (k[1] < 700 and k[1] > 5) ]

    max_coord = ccc
    
    
    for point in max_coord:
    
        startx = int(round(point[0] - bbox))
        endx = startx + bbox*2
        starty = int(round( point[1] - bbox))
        endy = starty + bbox*2

        if startx > 0 and startx < img.shape[0] - bbox*2:
            if starty > 0 and starty < img.shape[1] - bbox*2: 

                defim.resize(defim.shape[0] + 1, bbox*2, bbox*2)
                defim[defcount] = img[startx:endx, starty:endy]
                defcoord.resize(defcoord.shape[0] + 1, 2)
                defcoord[defcount] = point[0:2]

                defcount = defcount + 1
            
    
    return thresh, defim, defcoord

Here is an example of how it works on a single frame:

In [ ]:
thresh, defim, defcoord = extractDefects(imgdata[fnum], classpred, softmax_threhold = 0.95, bbox=16)
In [150]:
fig302 = plt.figure(302, figsize=(5,5))
ax = fig302.add_axes([0, 0, 1, 1])
ax.imshow(thresh)
bbox = 16

for point in defcoord:
    startx = int(round(point[0] - bbox))
    endx = startx + bbox*2
    starty = int(round( point[1] - bbox))
    endy = starty + bbox*2
    
    p = patches.Rectangle( (starty, startx), bbox*2, bbox*2,  fill=False, edgecolor='red', lw=4)
    ax.add_patch(p)
    
ax.axis('off')
Out[150]:
(-0.5, 511.5, 511.5, -0.5)

3.3 Applying to the whole movie

Now that we have the functions for the single frames, we can iterate over the movie and extract all the defects for classification

In [416]:
mdirectory = './Figures/' + filename[:-4] + '/Movie/'
if not os.path.exists(mdirectory):
    os.makedirs(mdirectory)
In [ ]:
bbox = 16

totdefim = np.zeros(shape=(0, bbox*2, bbox*2))
totdefcoord = np.zeros(shape=(0, 2))
totdeffnum = np.zeros(shape=(0,1))

fig303 = plt.figure(303, figsize=(5,5))
ax = fig303.add_axes([0, 0, 1, 1])


for i, img in enumerate(imgdata):
    #print(i)
    #1. use trained model to predict defecs
    classpred = predictDefects(img, model, target_size, nb_classes=2)
    
    #2. extract defects
    thresh, defim, defcoord = extractDefects(img, classpred, 0.95, bbox)
    
    #3. record the frame
    nd = len(defim)
    deffnum = np.zeros(shape=(nd,1)) + i
    
    #4. append to the array
    totdefim = np.append(totdefim, defim, 0)
    totdefcoord = np.append(totdefcoord, defcoord, 0)
    totdeffnum = np.append(totdeffnum, deffnum, 0)
    
    ax.cla()
    ax.imshow(img)
    for point in defcoord:
        startx = int(round(point[0] - bbox))
        endx = startx + bbox*2
        starty = int(round( point[1] - bbox))
        endy = starty + bbox*2

        p = patches.Rectangle( (starty, startx), bbox*2, bbox*2,  fill=False, edgecolor='red', lw=4)
        ax.add_patch(p)
    
    ax.axis('off')
    
    fig303.savefig(mdirectory + '_frame_' + str(i) + '.png')

plt.close(fig303)    
In [69]:
frames = [5, 25, 65, 85]
n = len(frames) + 1

plt.figure(figsize=(16, 8))
for i in range(1, n):
    ax = plt.subplot(2, n, i)
    i_ = frames[i-1]
    plt.imshow(imgdata[i_], cmap = 'gray')
    plt.title('Frame {0}\n'.format(i_) + 'Raw experimental data', fontsize = 10)
    ax.get_xaxis().set_visible(False)
    ax.get_yaxis().set_visible(False)
    plt.tight_layout()
    ax = plt.subplot(2, n, i + n)
    plt.imshow(imgdata[i_], cmap = 'gray')
    defcoord = totdefcoord[np.where(totdeffnum == i_)[0], :]
    for point in defcoord:
        startx = int(round(point[0] - bbox))
        endx = startx + bbox*2
        starty = int(round( point[1] - bbox))
        endy = starty + bbox*2
        p = patches.Rectangle( (starty, startx), bbox*2, bbox*2,  fill=False, edgecolor='red', lw=1)
        ax.add_patch(p)
    plt.title('Frame {0}\n'.format(i_) + 'Localized defects', fontsize = 10)
    ax.get_xaxis().set_visible(False)
    ax.get_yaxis().set_visible(False)
    plt.tight_layout()

Because the trained network relies only on the local characteristics of the image, it can identify defects on the later stages of system evolution when the long range periodicity of the lattice is broken due to second phase evolution and displacement and rotation of unreacted WS$_2$ fragments.

We can now save the data:

In [ ]:
np.save('defim.npy', totdefim)
np.save('defcoord.npy', totdefcoord)
np.save('deffnum.npy', totdeffnum)

4. Clustering the defects

In [11]:
totdefcoord = np.load('defcoord.npy')
totdeffnum = np.load('deffnum.npy')
totdefim = np.load('defim.npy')
bbox = 16

Now that we extracted the defects we can use gaussain mixture model method to cluster them

4.1 Estimating number of clusters

We need to find an estimation of number of different types of defects we can find in the movie. This is not a simple task, as we need to select number of clusters in such a way that it is both representative of most types we are interested in (e.g. not just bright vs dark), but, at the same time, does not start splitting hairs and overwhelming us with unnecessary details.

In [12]:
totdef = totdefim.reshape(totdefim.shape[0], totdefim.shape[1]*totdefim.shape[2])

4.1.1 Dendrogram

In [11]:
#taken from https://joernhees.de/blog/2015/08/26/scipy-hierarchical-clustering-and-dendrogram-tutorial/
def fancy_dendrogram(*args, **kwargs):
    max_d = kwargs.pop('max_d', None)
    if max_d and 'color_threshold' not in kwargs:
        kwargs['color_threshold'] = max_d
    annotate_above = kwargs.pop('annotate_above', 0)

    ddata = dendrogram(*args, **kwargs)
    
    if not kwargs.get('no_plot', False):
        plt.title('Hierarchical Clustering Dendrogram (truncated)')
        plt.xlabel('sample index or (cluster size)')
        plt.ylabel('distance')
        for i, d, c in zip(ddata['icoord'], ddata['dcoord'], ddata['color_list']):
            x = 0.5 * sum(i[1:3])
            y = d[1]
            if y > annotate_above:
                plt.plot(x, y, 'o', c=c)
                plt.annotate("%.3g" % y, (x, y), xytext=(0, -5),
                             textcoords='offset points',
                             va='top', ha='center')
        if max_d:
            plt.axhline(y=max_d, c='k')
    return ddata
In [12]:
Z = linkage(totdef , 'ward')
In [16]:
fig400 = plt.figure(400, dpi=100)
max_d = 33   # max_d as in max_distance
fancy_dendrogram(
    Z,
    truncate_mode='lastp',
    p=50,
    leaf_rotation=90.,
    leaf_font_size=6.,
    show_contracted=True,
    annotate_above=30,  # useful in small plots so annotations don't overlap
    max_d=max_d
)
plt.show()

Dendrogram is a one way to look at clustering. As distances between clusters become smaller, it makes less sense to plit them. In this example, we can estimate about 5-7 clusters with significant distances.

4.1.2 PCA

We can use 2 types of estimation from PCA. First, we can use dimensionality reduction method:

In [20]:
pca = PCA(n_components='mle', svd_solver = 'full')
In [21]:
pca.fit(totdef)
Out[21]:
PCA(copy=True, iterated_power='auto', n_components='mle', random_state=None,
  svd_solver='full', tol=0.0, whiten=False)
In [22]:
pca.n_components_
Out[22]:
250

In current run, it results in 250 components with significant contribution to variation. We can look at the Scree plot to find a more reasonable cut off:

In [24]:
plt.figure()
plt.plot(pca.explained_variance_ratio_[1:30], marker='o',markersize=4, linestyle='-', linewidth = 0.5, color = 'black')
plt.xlabel("number of components", fontsize = 12)
plt.ylabel('explained variance', fontsize = 12)
print(pca.explained_variance_ratio_[1:30]) 
[ 0.03939209  0.03510355  0.02098565  0.02052916  0.01766718  0.01363231
  0.01315034  0.01280795  0.01105607  0.00904456  0.0082742   0.00743088
  0.00701289  0.00669168  0.00574388  0.00515341  0.00507327  0.00472977
  0.00468524  0.00432358  0.00422907  0.00403801  0.00354088  0.00333545
  0.00302002  0.00288533  0.00282974  0.00276777  0.00264781]

In this case, plot begins to flatten out after about 5-15 components

4.1.3 BIC method

We can use Bayesian Information Criterion:

In [25]:
totdef2 = totdef - totdef.mean(axis=0)
In [26]:
maxc = 20
bics = list()
for i in range(maxc-1):
    gmm = mixture.GaussianMixture(n_components=i+2, covariance_type='full')
    gmm.fit(totdef2)
    bic = gmm.bic(totdef2)
    bics.append(bic)
    print(i +2, bic)
2 -10214721.2905
3 -6236494.16404
4 -2291641.83679
5 1570374.06464
6 5045864.16108
7 5757110.94101
8 8600335.14705
9 9015083.15442
10 10137010.8917
11 16035057.3676
12 17538349.8055
13 15411724.1257
14 17787512.0239
15 16905806.0022
16 25822572.7841
17 23571903.2442
18 27700862.4987
19 32528594.714
20 36488375.7737
In [27]:
plt.figure()
xt = np.arange(0,20,2)
plt.plot(bics, marker='o',markersize=4, linestyle='-', linewidth = 0.5, color = 'black')
plt.xlabel("number of components", fontsize = 12)
plt.ylabel("BIC", fontsize = 12)
plt.xticks(xt,fontsize = 12)
plt.yticks(fontsize = 12)
Out[27]:
(array([-20000000., -10000000.,         0.,  10000000.,  20000000.,
         30000000.,  40000000.]), <a list of 7 Text yticklabel objects>)

In this case, optimal number of clusters is 2. However, we can see some spikes in the graph at in between 5 and 15 clusters.

4.2 Identifying types of defects

Although the BIC results were somewhat ambiguous, we can choose the number of clusters we are interested in to be 5. We found it to be an optimal number of components for understanding the type of defect structures present in the data. Indeed, an increase in the number of components resulted in fine (sub-) structures of the detected defects, while decrease in the number of components produced some physically meaningless structures.

In [13]:
nc = 5

4.2.1 Clustering

We will use GMM to identify clusters and then construct an average image for each one:

In [14]:
gmm = mixture.GaussianMixture(n_components=nc, covariance_type='full', random_state=8000)
gmm.fit(totdef)
Out[14]:
GaussianMixture(covariance_type='full', init_params='kmeans', max_iter=100,
        means_init=None, n_components=5, n_init=1, precisions_init=None,
        random_state=8000, reg_covar=1e-06, tol=0.001, verbose=0,
        verbose_interval=10, warm_start=False, weights_init=None)
In [15]:
clusters = gmm.predict(totdef) + 1
In [16]:
colors = ['black', 'blue', 'red', 'cyan', 'magenta', 'green', 'darkviolet', 'gold', 'pink', 'red', 'yellow', 'maroon', 'orangered', 'lime', 'darkblue']
In [17]:
cla = np.ndarray(shape=(np.amax(clusters), 32,32))
fig401 = plt.figure(401, figsize=(10, 5))

rows = int(np.ceil(float(nc)/5))
cols = int(np.ceil(float(np.amax(clusters))/rows))

gs1 = gridspec.GridSpec(rows, cols)

for i in range(np.amax(clusters)):
    cl = totdefim[clusters == i + 1]
    cla[i] = np.mean(cl, axis=0)
    
    ax = fig401.add_subplot(gs1[i])
    ax.imshow(cla[i], cmap = 'gray', Interpolation = 'Gaussian')
    
    ax.axis('off')
    ax.set_title('Class ' + str(i + 1) + '\nCount: ' + str(len(cl)))
    p = patches.Rectangle( (0, 0), bbox*2-1, bbox*2-1,  fill=False, color = colors[i], lw=4)
    ax.add_patch(p)

#plt.subplots_adjust(hspace=0.6, wspace=0.4)
fig401.savefig(directory + 'gmm_nc_' + str(nc) + '.png')

4.2.2 Show clustering on the movie

In [18]:
mdirectory2 = './Figures/' + filename[:-4] + '/MovieClustered_reduced/'
if not os.path.exists(mdirectory2):
    os.makedirs(mdirectory2)
In [ ]:
fig403 = plt.figure(403, figsize=(5,5))
ax = fig403.add_axes([0, 0, 1, 1])


for i, img in enumerate(imgdata):
    
    ax.cla()
    ax.imshow(img, cmap ='gray')
    
    
    defcoord = totdefcoord[np.where(totdeffnum == i)[0], :]
    clust = clusters2[np.where(totdeffnum == i)[0]]
    
    for i2, point in enumerate(defcoord):
        startx = int(round(point[0] - bbox))
        endx = startx + bbox*2
        starty = int(round( point[1] - bbox))
        endy = starty + bbox*2

        
        
        p = patches.Rectangle( (starty, startx), bbox*2, bbox*2,  fill=False, edgecolor=colors[clust[i2] - 1], lw=4)
        ax.add_patch(p)
    
    ax.axis('off')
    fig403.savefig(mdirectory2 + 'frame_' + str(i) + '.png')

plt.close(fig403)

Finally, we should save the data:

In [370]:
#save the clustering data
np.save('clusters.npy', clusters)
np.save('clusters_reduced.npy', clusters2)

4.3 Distribution of defects in time

4.3.1 Total counts

We are going to calculate total number of defects and their class distribution in each frame

To view the trends better, we can use moving average to smooth the curves:

In [19]:
def moving_average(a, n=3) :
    ret = np.cumsum(a, dtype=float)
    ret[n:] = ret[n:] - ret[:-n]
    return ret[n - 1:] / n
In [20]:
nframes = imgdata.shape[0]
totclust = np.zeros(shape = (nframes))
totbyclust = np.zeros(shape = (nc, nframes))
totbyclust2 = np.zeros(shape = (nc, nframes))
for i in range(nframes):
    clust = clusters[np.where(totdeffnum == i)[0]]
    totclust[i] = len(clust)
    
    for i2 in range(nc):
        totbyclust[i2, i] = (clust == i2 + 1).sum()    
In [30]:
fig430 = plt.figure(430, dpi=150)
ax = fig430.add_axes([0.1, 0.1, 0.7, 0.9])
ax.plot(moving_average(totclust, 5))
ax.set_xlim(0, 99)
ax.set_xlabel('Frame', fontsize = 14)
ax.set_ylabel('Defect Count', fontsize = 14)
ax.tick_params(axis='both', which='major', labelsize=11)

Here we can observe a significant increase in defects starting at around frame 50. If we look a the movie, this is about when the holes start appearing. Total number of defects peaks at frame 80, due to holes in the material occupying about 25% of the image.

4.3.2. "Broom" graph

To better understand time evolution of defects, we can construct graphs tracking defects types in (x, y, t) coordinate system

In [32]:
def BroomLines(spo):
    lcin = np.ndarray(shape=(len(spo), 2, 3))
        
    lcin[:, 0, 2] = spo[:, 0]
    lcin[:, 1, 2] = spo[:, 0]
    lcin[:, 0, 1] = spo[:, 1]
    lcin[:, 1, 1] = spo[:, 1]
    lcin[:, 0, 0] = i
    lcin[:, 1, 0] = (i+1) 
    
    return lcin
    
In [33]:
fig432 = plt.figure(432, figsize=(8, 8))
ax = fig432.add_subplot(111, projection='3d')
# We will look into behavior of defects associated with 3 distorted sattes of Mo_w dopant
nc_ = [1, 2, 3, 4, 5]

for i in range(nframes):
    clust = clusters[np.where(totdeffnum == i)[0]]
    coord = totdefcoord[np.where(totdeffnum == i)[0]]

    for c, i2 in enumerate(nc_):
        spo = coord[clust==i2]
        ctest2 = BroomLines(spo)
        lc = Line3DCollection(ctest2, linewidths=1.5, colors=colors[c])    
        ax.add_collection3d(lc)
        
        
ax.set_zlim3d([0, 512])
ax.axis([0, 100, 0, 512])
        
ax.set_zlabel('X')
ax.set_ylabel('Y')
ax.set_xlabel('Frame')
Out[33]:
<matplotlib.text.Text at 0x22f57fde080>

There are three statistical behaviors: weakly moving trajectories, stronger diffusion, and “uncorrelated events”. The most well-defined trajectories are associated with Mo dopants (class 1 and class 3). Interestingly, these Mo defects show different diffusion behaviors depending on their location in the lattice and are characterized by reversible switching between two configurations (class 1 and class 3) along their trajectories. The difference in diffusion behavior may be potentially connected to complex spatial character of strain fields during the material transformation, which may impact diffusion properties

We can further collapse the representation to 2D. In order to show differences in time, we can project windows of defects we found. This allows us to show defects that are continous in time more clearly than randomly occuring ones.

In [22]:
imsum = np.zeros(shape = (imgsrc.shape[0],imgsrc.shape[1]) )

fig435 = plt.figure(435, figsize=(8, 8))
ax = fig435.add_subplot(111)

#list all the clusters for plotting 
cltoplot = [5]

coord = np.empty((0, 2))
for cl in cltoplot:
    coord_ = totdefcoord[np.where(clusters == cl)[0]]
    coord = np.concatenate((coord, coord_))
    for i2, point in enumerate(coord):
        startx = int(round(point[0] - bbox))
        endx = startx + bbox*2
        starty = int(round( point[1] - bbox))
        endy = starty + bbox*2

        imsum[starty:endy, startx:endx] = imsum[starty:endy, startx:endx] + 1

immean = imsum/imgdata.shape[0]

imsh = ax.imshow(immean, cmap='nipy_spectral', Interpolation = 'Gaussian')
plt.colorbar(imsh)
plt.ylim(0, 512)
Out[22]:
(0, 512)

We will now isolate selected trajectories for future analysis of diffusion and transition between different states

In [23]:
def rem_coord(co, center, dist_center):
    return [co[0] < center[0] - dist_center, co[0] > center[0] + dist_center,
            co[1] < center[1] - dist_center, co[1] > center[1] + dist_center] 

nc_ = [5]
#selecting projections with highest intensity in the analysis based on projected windows
rois = [(256, 283), (39, 123)]
flow = {}
flow_2d = {}
for i_c, center in enumerate(rois):
    cl_coord = {}
    cl_coord_all_fr = np.empty((0, 2))
    for i in range(nframes):
        clust = clusters[np.where(totdeffnum == i)[0]]
        coord = totdefcoord[np.where(totdeffnum == i)[0]]
        cl_coord_fr = OrderedDict()
        for i2 in nc_:
            spo = coord[clust==i2]
            coord_to_rem = [idx for idx, c in enumerate(spo) if any(rem_coord(c, center, 25))]
            coord_to_rem = np.array(coord_to_rem, dtype = int)
            spo = np.delete(spo, coord_to_rem, axis = 0)
            cl_coord_fr[i2] = spo
            cl_coord_all_fr = np.concatenate((cl_coord_all_fr, spo))          
        cl_coord[i] = cl_coord_fr
    flow_2d[i_c] = cl_coord_all_fr
    flow[i_c] = cl_coord
In [24]:
#np.save("flow_S_cl5.npy", flow)

Diffusion and switching behaviour for selected trajectories

We will need these functions

In [25]:
def diff_coef(x, y):
    ra = np.var(x)
    rb = np.var(y)
    r = ra + rb
    D = r/(4*len(x)*frame_t)
    return D
    
def eigsorted(cov):
    vals, vecs = np.linalg.eigh(cov)
    order = vals.argsort()[::-1]
    return vals[order], vecs[:,order]

We can now plot 2-d projections of selected "flows" with calculated prediction ellipses. We will also estimate a diffusion coefficient usingn a simple 2-d random walk model.

In [26]:
# confidence interval
nstd = 2
# conversion coefficient (pixels to nanometers) 
px_to_nm = 0.03
# conversion coefficient (frames to time(s))
frame_t = 2.25
In [28]:
from matplotlib.patches import Ellipse
In [29]:
fig436 = plt.figure(436, figsize = (5, 5), dpi = 150)
ax = fig436.add_subplot(111, aspect='equal')
ax.set_xlabel('X', style = 'italic', fontsize = 14)
ax.set_ylabel('Y', style = 'italic', fontsize = 14)
for k_f in range(len(flow_2d.keys())):
    x = flow_2d[k_f][:, 0]*px_to_nm
    y = flow_2d[k_f][:, 1]*px_to_nm
    plt.plot(x, y, 'o', c = 'magenta', markersize = 14, alpha = 0.6, zorder = 1)
    cov = np.cov(x, y)
    vals, vecs = eigsorted(cov)
    theta = np.degrees(np.arctan2(*vecs[:,0][::-1]))
    w, h = 2 * nstd * np.sqrt(vals)
    ell = Ellipse(xy=(np.mean(x), np.mean(y)),
                  width=w, height=h,
                  angle=theta, color='black', linewidth = 1.25, ls = 'dotted', zorder = 2)
    ell.set_facecolor('none')
    ax.add_artist(ell)
    ax.text(np.mean(x), np.mean(y), str(k_f) , fontsize=6)
    print('Diffusion coefficient for trajectory ' + str(k_f), ':', diff_coef(x, y), 'nm^2/s')    
Diffusion coefficient for trajectory 0 : 0.000289121441816 nm^2/s
Diffusion coefficient for trajectory 1 : 0.000625766184957 nm^2/s

We can also make a 1-d r(t) representation of the selcted "flows"

In [30]:
fig = plt.figure(figsize = (5, 5), dpi = 150)
ax = fig.add_subplot(111)
ax.set_xlabel('Time (s)', fontsize = 16)
ax.set_ylabel('r', style = 'italic', fontsize = 16)
ax.tick_params(axis='x', labelsize=10)
ax.tick_params(axis='y', labelsize=10)
n_frames = 100
for k_f in range(len(flow.keys())):
    spo = flow[k_f]
    for frame in range(n_frames):
        for idx, (k, v) in enumerate(spo[frame].items()):
            if v.shape[0] != 0: 
                v_r = np.sqrt((v[0][0]*px_to_nm)**2 + (v[0][1]*px_to_nm)**2)
                plt.plot(frame*frame_t, v_r, 'o', c=colors[k-1], markersize = 8, alpha = 0.6)  
In [ ]: